# Project Info---------------
#  Project: Outgroup Avoidance
#  Purpose: Analyze Second Survey
#  Outputs: Figure 9, S7, Table S6-7 and S10

#Load packages and data-----

library("tidyverse")
library("estimatr")
library("stdidx")
library("modelsummary")
library("knitr")
library("kableExtra")

stdy3_dta <- read.csv("study3.csv")

#Main Text------
##Figure 9-----
###Panel A-----
#prop.table(table(stdy3_dta$y_choose_pal[stdy3_dta$d_any_treatment==0]))#0.14 of sample select pal posts
stdy3_dta %>% 
  filter(.,
         !is.na(y_choose_pal),
         d_any_treatment ==0) %>% 
  ggplot(., aes(x=factor(y_choose_pal)))+
  geom_bar(width=0.7, fill="dodgerblue4")+
  scale_x_discrete(labels = c("No", "Yes"))+
  xlab("Select Palestinian Post")+
  ylab("Count")+
  theme_bw()+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 18),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())

###Panel B-----
pooled <-
  lm_robust(y_choose_pal ~ d_any_treatment,
            data = stdy3_dta) %>% 
  tidy(.) %>% 
  filter(.,
         term != "(Intercept)") %>% 
  mutate(.,
         model = "Pooled",
         term = case_when(
           term == "d_any_treatment" ~ "All Treatments"))
#order treatments so that the control group is a reference category
stdy3_dta$treat <-
  relevel(as.factor(stdy3_dta$treat) , ref = "control")

by_treatment <-
  lm_robust(y_choose_pal ~ treat,
            data = stdy3_dta)%>% 
  tidy(.) %>% 
  filter(.,
         term != "(Intercept)") %>% 
  mutate(.,
         model = "By Treatment",
         term = case_when(
           term == "treataffirmation" ~ "Affirmation",
           term == "treatendorsement" ~ "Ingroup Endorsment",
           term == "treatfact" ~ "Fact Checking",
           term == "treathumility" ~ "Intellectual Humility",
           term == "treatincentive" ~ "Monetary Incentive",
           term == "treatmalleability" ~ "Malleability",
           term == "treatplacebo" ~ "Placebo"
         ),
         p.value = p.adjust(p.value, method = "BH"))


#Plot Estiamts
main_results <- 
  rbind(pooled, by_treatment) %>% 
  mutate(.,
         b = paste0("b=", round(estimate, 3)),
         p = paste0("p=", round(p.adjust(p.value, method = "BH"), 3)),
         bp = paste0(b, ", ", p))

main_results$model <- 
  factor(main_results$model, 
         c("Pooled", "By Treatment"))



ggplot(main_results, aes(x=estimate, y=term)) + 
  geom_point(position = position_dodge(width = 0.6))+
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high),
                position = position_dodge(width = 0.6),
                width=.1) +
  ggforce::facet_col(facets = vars(model), 
                     scales = "free_y", 
                     space = "free") +
  geom_vline(xintercept = 0,
             color = "black",
             linetype = "dashed",
             size = 0.4) +
  #        geom_label()+
  geom_text(aes(label= bp), 
            colour = "dodgerblue4",
            position = position_dodge(width = .6),
            hjust = .4,
            vjust = -0.6,
            size = 4,
            show.legend = FALSE) +
  labs(x = "",
       y = "",
       color = "")+
  #  xlim(-.25, .5)+
  scale_color_manual(values = c("dodgerblue3"))+
  # theme_tufte()+
  guides(color = guide_legend(reverse = TRUE))+
  theme_bw()+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 18),
        # axis.text.y = element_text(size=14),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())



#Appendix----
##Figure S7----
stdy3_dta$placebo_ref <-
  relevel(as.factor(stdy3_dta$treat) , ref = "placebo")

stdy3_dta %>% 
  filter(.,
         placebo_ref!="control") %>% 
  lm_robust(y_choose_pal ~ placebo_ref,
            data = .) %>% 
  tidy(.) %>% 
  filter(.,
         term != "(Intercept)") %>% 
  mutate(.,
         model = "By Treatment",
         term = case_when(
           term == "placebo_refaffirmation" ~ "Affirmation",
           term == "placebo_refendorsement" ~ "Ingroup Endorsment",
           term == "placebo_reffact" ~ "Fact Checking",
           term == "placebo_refhumility" ~ "Intellectual Humility",
           term == "placebo_refincentive" ~ "Monetary Incentive",
           term == "placebo_refmalleability" ~ "Malleability",
           term == "placebo_refpersonal" ~ "Personal Content",
           term == "placebo_refplacebo" ~ "Placebo"
         ),
         p.value = p.adjust(p.value, method = "BH")) %>% 
  mutate(.,
         b = paste0("b=", round(estimate, 3)),
         p = paste0("p=", round(p.value, 3)),
         bp = paste0(b, ", ", p)) %>% 
ggplot(., aes(x=estimate, y=term)) + 
  geom_point(position = position_dodge(width = 0.6))+
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high),
                position = position_dodge(width = 0.6),
                width=.1) +
  geom_vline(xintercept = 0,
             color = "black",
             linetype = "dashed",
             size = 0.4) +
  #        geom_label()+
  geom_text(aes(label= bp), 
            colour = "dodgerblue4",
            position = position_dodge(width = .6),
            hjust = .4,
            vjust = -0.5,
            size = 4,
            show.legend = FALSE) +
  labs(x = "",
       y = "",
       color = "")+
  scale_color_manual(values = c("dodgerblue3"))+
  guides(color = guide_legend(reverse = TRUE))+
  theme_bw()+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 18),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())

##Table S6-----
affirm <- 
  lm_robust(mnp_affirm ~ d_affirm,
            data = stdy3_dta)

humil <-
  lm_robust(mnp_humil ~ d_humil,
            data = stdy3_dta)

mall <-
  lm_robust(mnp_mall ~ d_mall,
            data = stdy3_dta)

means <- 
  data.frame(t(c("Control Mean",
                 round(mean(stdy3_dta$mnp_affirm[stdy3_dta$treat=="control"], na.rm = T), 3),
                 round(mean(stdy3_dta$mnp_humil[stdy3_dta$treat=="control"], na.rm = T), 3),
                 round(mean(stdy3_dta$mnp_mall[stdy3_dta$treat=="control"], na.rm = T), 3))))

sds <- 
  data.frame(t(c("Control SD",
                 round(sd(stdy3_dta$mnp_affirm[stdy3_dta$treat=="control"], na.rm = T), 3),
                 round(sd(stdy3_dta$mnp_humil[stdy3_dta$treat=="control"], na.rm = T), 3),
                 round(sd(stdy3_dta$mnp_mall[stdy3_dta$treat=="control"], na.rm = T), 3))))

meansd <- rbind(means,sds)

models <- 
  list("Affirmation Index" = affirm, 
       "Humility Index" = humil, 
       "Malleability Index" = mall)
sig_start <- c('*' = .05) 

modelsummary(models,
             output ="latex",
               title = "Average Treatment Effect on Psychological Constructs\\label{manip_study3}",
               coef_map = c(d_affirm = "Affirmation Treatment", d_humil = "Int. Humility Treatment",
                            d_mall = "Malleability Treatment"),
               gof_map = c("nobs", "r.squared"),
               stars = sig_start,
               add_rows = meansd)  




##Table S7-----
tot_mall <-
  iv_robust(y_choose_pal ~ 
              mnp_mall|
              d_mall,
            data = stdy3_dta)

tot_aff <-
  iv_robust(y_choose_pal ~ 
              mnp_affirm |
              d_affirm,
            data = stdy3_dta)

tot_humil<-
  iv_robust(y_choose_pal ~ 
              mnp_humil |
              d_humil,
            data = stdy3_dta)


outcome_means <- 
  data.frame(t(c("Control Mean",
                 round(mean(stdy3_dta$y_choose_pal[stdy3_dta$treat=="control"], na.rm = T), 3),
                 round(mean(stdy3_dta$y_choose_pal[stdy3_dta$treat=="control"], na.rm = T), 3),
                 round(mean(stdy3_dta$y_choose_pal[stdy3_dta$treat=="control"], na.rm = T), 3))))

outcome_sds <- 
  data.frame(t(c("Control SD",
                 round(sd(stdy3_dta$y_choose_pal[stdy3_dta$treat=="control"], na.rm = T), 3),
                 round(sd(stdy3_dta$y_choose_pal[stdy3_dta$treat=="control"], na.rm = T), 3),
                 round(sd(stdy3_dta$y_choose_pal[stdy3_dta$treat=="control"], na.rm = T), 3))))

outcome_meansd <- rbind(outcome_means,outcome_sds)

tot_models <- 
  list("1" = tot_aff, 
       "2" = tot_humil, 
       "3" = tot_mall)
sig_start <- c('*' = .05) 

modelsummary(tot_models,
               output ="latex",
               title = "TOT Estimates for Study 3 Psychological Treatments\\label{tot_study3}",
               coef_map = c(mnp_affirm = "Affirmation Treatment", mnp_humil = "Int. Humility Treatment",
                            mnp_mall = "Malleability Treatment"),
               gof_map = c("nobs", "r.squared"),
               stars = sig_start,
               add_rows = outcome_meansd) %>%
  add_header_above(c("Outcome: Select Palestinian Post (0/1)" = 4))


##Table S10-----
covs <- 
  stdy3_dta %>%
  mutate(.,
         Age = case_when(
           age == "1" ~ "18-24",
           age == "2" ~ "25-34",
           age == "3" ~ "35-44",
           age == "4" ~ "45-54",
           age == "5" ~ "55-64",
           age == "6" ~ "65-70"
         ),
         Gender = case_when(
           gnd == "1" ~ "Male",
           gnd == "2" ~ "Female"
         ),
         Religiosity = case_when(
           rlg == "1" ~ "Secular",
           rlg == "2" ~ "Traditional",
           rlg == "3" ~ "Religious",
           rlg == "4" ~ "Ultra-Orthodox"
         )) %>% 
  dplyr::select(Age, Gender, Religiosity)

# Calculate proportions excluding "NAs"
proportions_df <- 
  covs %>%
  pivot_longer(cols = everything(), names_to = "Variable", values_to = "Value") %>%
  filter(!is.na(Value)) %>%  # Exclude "No-Response"
  group_by(Variable, Value) %>%
  summarise(
    Count = n(),  # Count of each value
    .groups = "drop_last"  # Retain grouping by Variable for next step
  ) %>%
  mutate(
    Total = sum(Count),  # Total counts for each variable
    Proportion = Count / Total  # Proportion for each value
  ) %>%
  ungroup() %>%  # Remove all grouping
  arrange(Variable, desc(Proportion)) %>% 
  dplyr:: select(Variable,Value,Proportion) %>% 
  mutate(.,
         Proportion = round(Proportion,3))


# read in ipanel data
ipanel <- read.csv("ipanel_quota.csv")
`Ipanel prop.` <-ipanel$IPanel.Proportion
descriptive_table <- cbind(proportions_df,`Ipanel prop.`)

# Save the LaTeX table
descriptive_table %>%
  kable(format = "latex", 
        booktabs = TRUE, 
        align = "c", 
        label = "tab:desc_3",
        caption = "Descriptive Statistics Study 3") 
